clear all
close all
clc

iteration=1;
ae=1;
I=rgb2gray(imread('IMG_0002.jpg'));




I3=double(I>200);

I32=double(I>200);
% I=imread('216.jpg');
% I=bwmorph(I,'open');
% I=bwmorph(I,'close');
% I3=double(I);
% I32=double(I);
test=0;
 while iteration<4
CC=bwconncomp(I3);
numPixels = cellfun(@numel,CC.PixelIdxList);
[biggest,idx] = max(numPixels); 
% I2=I;
% I2(CC.PixelIdxList{idx}) = 2;
% I2=I2>1;
% l=CC.PixelIdxList{idx};
% [x y]=ind2sub(size(I),l);
[ax ay]=ind2sub(size(I3),cell2mat(CC.PixelIdxList(idx)));
coord=[ax ay];
I3(cell2mat(CC.PixelIdxList(idx)))=2;
I3=I3>1;
P=regionprops(CC,'all');
C=P(idx).Centroid;
if iteration==1;
    compactness=P(idx).Area/(4*pi*(P(idx).Perimeter)^2);
    if compactness<0.0016
        result='dear'
    else if compactness <0.0021
            result='squirell'
        else if compactness <0.0027
                result='fox'
            else
                result='turkey'
            end
        end
    end
end
            

ed=edge(I3);
CC2=bwconncomp(ed);
azer=[];
 for k=1:CC2.NumObjects
azer=[azer; cell2mat(CC2.PixelIdxList(k))];
end
[ax2 ay2]=ind2sub(size(ed),azer);
edx=[ax2 ay2];
distanc=pdist2(coord,edx);
for k = 1:size(distanc,1)
% dd(k,:)=sort(distanc(k,:));
ss(k)=sum(diff(sort(distanc(k,:))));
sss(k)=sum(sort(distanc(k,:)).*[ones(1,20) zeros(1,size(distanc,2)-20)]);

%  ssss(k)=sum(dd(k,1:20));
end
% ss=median(distanc,2);
[ccccc indic]=min(ss./(sss));
C(1)=coord(indic,2);
C(2)=coord(indic,1);
D=[];
% i=0;
% skel=bwmorph(I,'skel',Inf);
% I2=I>200;

            for k=1:size(ed,1)
                for l=1:size(ed,2)
                    if ed(k,l)
                        d=sqrt((C(1)-l)*(C(1)-l)+(C(2)-k)*(C(2)-k));
                        d2=[d k l];
                        D=[D ;d2];
                    end
                end
            end
            m=min(D(:,1));
            ind=find(D(:,1)==min(D(:,1)),1);
            ymin=D(ind,2);
            xmin=D(ind,3);           






v=[xmin-C(1) ymin-C(2)];
v2=[-v(2) v(1)]/sqrt(v(1)*v(1)+v(2)*v(2));
v2=v2/sqrt(v2(1)*v2(1)+v2(2)*v2(2));
a=0;

e=0;
e2=0;

while (e~=1&e2~=1)
    a=1+a;
%     if (floor(C(2)+a*v2(2))<1|floor(C(1)+a*v2(1))<1)
%         e=1;
%     else if (floor(C(1)-a*v2(1))<1|floor(C(2)-a*v2(2))<1)
%             e2=1;
%         else
    e=ed(floor(C(2)+a*v2(2)),floor(C(1)+a*v2(1)));
    e2=ed(floor(C(2)-a*v2(2)),floor(C(1)-a*v2(1)));
%         end
%     end
end

if e==1
xmax=C(1)+(a)*v2(1);
ymax=C(2)+(a)*v2(2);
else

xmax=C(1)-(a)*v2(1);
ymax=C(2)-(a)*v2(2);
end
major=sqrt((C(1)-xmax)*(C(1)-xmax)+(C(2)-ymax)*(C(2)-ymax));

% ang=-(P(1).Orientation+90)*pi/180;
if (ymax-C(2))>0
ang=acos((xmax-C(1))/major)-pi/2;
else
    ang=pi/2-acos((xmax-C(1))/major);
end


t=pi/2-acos((xmax-C(1))/major);
if (ymax-C(2))<0

    t=acos((xmax-C(1))/major)-pi/2;
end
 

a=major;
b=min(D(:,1));

for x=1:size(I3,1)
    for y=1:size(I3,2)
 X = (x-C(2))*cos(t)+(y-C(1))*sin(t);% Translate and rotate coords.
 Y = -(x-C(2))*sin(t)+(y-C(1))*cos(t); % to align with ellipse
if ((X*X)/(a*a)+(Y*Y)/(b*b))<1
    I32(x,y)=0;
    
end
    end
end

% iteration=sum(sum(abs(I32-I3)));
iteration=iteration+1;
I3=I32;
ae=ae+1;
test=test+1;
% if ae==1
figure()
% subplot(1,3,iteration-1)
imshow(I3,[])
% end
hold on
plot(C(1),C(2),'ob');
plot([C(1) xmin],[C(2) ymin],'r');
plot([C(1) xmax],[C(2) ymax],'y');
h=ellipse(min(D(:,1)),min(D(:,1)),ang,C(1),C(2),'r');
  g=ellipse(min(D(:,1)),major,ang,C(1),C(2),'y');
clear distanc ss sss
 end